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This article presents a new control approach and a dynamic model for engineered flap- 
ping flight with many interacting degrees of freedom. This paper explores the applications 
of neurobiologically inspired control systems in the form of Central Pattern Generators 
(CPG) to control flapping flight dynamics. A rigorous mathematical and control theoretic 
framework to design complex three dimensional wing motions is presented based on phase 
synchronization of nonlinear oscillators. In particular, we show the flapping flying dynamics 
without a tail or traditional aerodynamic control surfaces can be eff*ectively controlled by 
a reduced set of CPG parameters that generate phase-synchronized or symmetry-breaking 
oscillatory motions of two main wings. Furthermore, by using Hopf bifurcation, we show 
that tailless aircraft alternating between flapping and gliding can be effectively stabilized 
by smooth wing motions driven by the CPG network. Results of numerical simulation with 
a full six degree-of-freedom flight dynamic model validate the effectiveness of the proposed 
neurobiologically inspired control approach. 

Nomenclature 



4>w^ i^w, Ow Flapping, lead-lag, and pitch angles of each wing (left, right) 

= {ui,Vi)'^ State vector of the i-th Hopf oscillator 

f (x^; pi) Hopf nonlinear equations in the vector form with radus pi 

Pi Radius of the limit cycle from the i-th Hopf oscillator 

A Common rate of convergence of Hopf oscillators 

uj Common oscillation frequency of Hopf oscillators, rad/s 

tti Amplitude bias of the i-th Hopf oscillator 

(7 Bifurcation parameter, a = 1 for a stable limit cycle or cr = —1 for convergence to a^. 

R(Aij) 2x2 rotational transformation matrix 

Aij Phase lead of the i-th Hopf oscillator from the j-th 

n Total number of Hopf oscillators in the CPG network 

Ijt Identity matrix G M^^^ 

k Coupling gain of the coupled Hopf oscillators 

kr Reduced frequency of the flapping wing 

£ Contraction rate of the virtual nonlinear system 

G Original Laplacian matrix with rotational transformation G M?^ ^ -^^ 

L Graph Laplacian matrix G M2nx2n 

V Matrix of orthonormal eigenvectors of L without the ones vector G '^'^'^^'^i'^—'^) 

^min(-), Amacc(-) Minimum or maximum eigenvalues of the matrix 

= (xb^Vb^ ^b) Vehicle body frame coordinates 

y^w = {xw,yw: Zs) Wing frame coordinates (left, right) 

Xs = (xs,ys,Zs) Stroke plane frame coordinates 

R Wing span of a single wing, m 

r Wing span coordinate value r G [0, i?] of the wing blade element, m 

Ciwir, t) Local angle of attack of the wing blade element, rad 

/5it;(T, t) Local direction of the wind in the wing frame, rad 

Vfe Speed of the vehicle, without the relative wind, m/s 

Vr Local wind speed of the wing blade element, m/s 



= {Vwx^ Vwy, Vwz)^ Total wind velocity vector of the blade element, m/s 
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Local lift and drag coefficients of the blade element 

Total aerodynamic forces of each wing in the x and z directions of the wing frame 
Total aerodynamic forces from the wing in the body frame (left, right) 
Location of the stroke and wing frames with respect to the body frame 
Angular transformation matrix 

Angle of attack and slide-slip angle of the vehicle body, rad/s 
Stroke plane inclination angle from the vertical line 
Directional cosine matrix of Euler angles G ^ ^ 
Body angular rates, rad/s 

Euler angles of the vehicle body with respect to the inertial frame, rad 
Mass of the aircraft, kg 

Total aerodynamic moments from each wing, Nm 

Vehicle velocity vector in the body frame, m/s 

Gravitational force vector in the inertial frame 

Additional forces generated by the body (fuselage) and the tail 

Inertia matrix G M"^^"^ , kgm^ 

Additional aerodynamic moments from the body and the tail 




m 



m={Mx,My,Mz)^ 

Vfe = (y,,,y,^,H,)^ 

= (0,0, mp)^ 

h 

B = {Bx,By,Bz)^ 

Subscript 



right, left 
6, s, w 



Variable number of the coupled Hopf nonlinear oscillators 

Right or left wing 

Body, stroke, and wing frame 



I. 



Introduction 



Engineered flapping flight holds promise for creating biomimetic micro aerial vehicles (MAVs) flying in low 
Reynolds number regimes (Re< 10^) where rigid fixed wings drop substantially in aerodynamic performance. 
MAVs are typically classified as having maximum dimensions of 15 cm and fiying at a nominal speed of 1-20 
m/s in tight urban environments.^'^ Although natural flyers such as bats, birds, and insects have captured 
the imaginations of scientists and engineers for centuries, the maneuvering characteristics of unmanned aerial 
vehicles (UAVs) are nowhere near the agility and efficiency of animal ffight.^"^ Such highly maneuverable 
MAVs will make paradigm-shifting advances in monitoring of critical infrastructures such as power grids, 
bridges, and borders, as well as in intelligence, surveillance, and reconnaissance applications. 

The objective of this article is to investigate whether the control and synchronization of coupled non- 
linear oscillators, inspired by central pattern generators (CPGs) found in animal spinal cords, can control 
biomimetic flapping flight (see Fig. 1). An engineered CPG network, which ensures the stability and robust 
adaptation of motion, can significantly reduce the complexity associated with fiapping fiight. Unique to this 
research approach is the potential to reverse-engineer the key mechanisms of highly adaptive and robust 
rhythmic pattern modulations of flapping flight by integrating the neurobiological principles with the rigor- 
ous mathematical tools borrowed from nonlinear synchronization theory^ and fiight dynamics and controls. 
Such an approach has not been adopted for engineered flapping flight. 

While unsteady aerodynamics of flapping flight in low Reynolds number regimes has been extensively 
studied through numerical^' ''"^^ and experimental studies,^' one of the most interesting and least un- 
derstood aspects of flapping flight is how to precisely control and synchronize a large number of interacting 
limbs and joints that generate complex three-dimensional oscillatory movements of the wings governed by 
unsteady aerodynamic forces. There is relatively less prior work in control of fiapping fiight with notable 
exceptions of [14,18-21]. In this paper, we focus on three stereotyped motion primitives to define the three 
dimensional movements of wings: main fiapping (stroke) motion (Fig. 2a), horizontal wing-sweeping mo- 
tion called lead-lag motion (Fig. 2b), and wing pitch twisting (Fig. 2c). Studying how to produce such 
synchronized wing motions is expected to shed light on the key characteristics of animal fiapping fiyers. 

Previous robotic fiapping fiyers and their control design consider one or two degrees of freedom in the 
wings.-^^'^^'-^^'-^^'^^'^^"^^ However, even insects like the dragonfiy {Anax parthenope) are reported to have 
complex three-dimensional movements by actively controlling fiapping and twisting of four independent 
wings. ^ 

Furthermore, prior studies in fiapping flight^"^'^'^"^^'^^'^'"'^^"^^ assumed a very simple sinusoidal function 
for each joint to generate fiapping oscillations, without deliberating on how multiple limbs (or their nervous 
systems) are connected and actuated to follow such a time- varying reference trajectory. However, as shall 
be seen later in this paper, the use of sinusoidal functions to generate oscillatory motions of the wings is 
not effective for stable and agile ffapping ffying maneuvers especially with time- varying oscillation frequency. 
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Figure 1. Proposed hierarchical control structure with a CPG network with a reduced set of commands such as 
frequency uj{t), phase difference between wing joint angles A^^ , and max amplitude pi. 




(a) flapping 



(b) lead-lag sweeping {iJjw) 



(c) pitch {6w) and cambering 



Figure 2. Basic wing movements of bats employed in this paper (pictures from [17]). Except for cambering, birds 
exhibit similar wing movements. 



amplitude, and bias. Experimental results using high speed cameras have shown that the flapping motions in 
bats and birds are more complicated than perfect sinusoidal^' with a flxed amplitude. Furthermore, we show 
in this paper that phase differences between multiple oscillators, which are not permitted by independent 
control of each wing joint, can be effective control mechanisms for flapping flight dynamics. In order to 
bridge this gap, this article aims to establish a novel adaptive CPG-based control theory for flapping flight, 
through neuromechanical modeling, nonlinear control and synchronization, and numerical simulation. 

The paper is organized as follows. We illustrate the fundamentals and advantages of the CPG based con- 
trol for engineered flapping flight in Section II. We present a mathematical and control theoretic formulation 
of synchronized motions of multiple joints in the wings and body in Section II. C. The kinematic modeling of 
three-dimensional multi-joint wings and the dynamic models of flapping flight are presented in Section III. 
In Section IV, we present key control strategies and results of simulation with multijoint synchronization 
that enables a smooth transition between stable flapping flight and gliding flight. Concluding remarks are 
presented in Section V. We understand the challenges associated with building lightweight actuators that 
must be overcome to truly realize the potential of three dimensional wing maneuvers. We present the funda- 
mental neurobiologically inspired control theory that can further contribute to engineered flapping flight once 
such light-weight actuators become available in the future. In the meantime, we show how the multi-joint 
robotic bat testbed (see Fig. 4) driven by CPG control can further enhance our understanding of biomimetic 
flapping flight. 



II. Fundamentals of Neurobiologically Inspired Control 

This article reports the flrst investigation of CPG models by using coupled limit cycle oscillators for the 
purpose of controlled engineered flapping flight. Also, we introduce the use of Hopf bifurcation to enable 
smooth switching between flapping and gliding flight. Hooper^^ deflnes the central pattern generators of 
animals as neural networks that can endogenously (i.e., without rhythmic sensory or central brain input) 
produce coordinated patterns of rhythmic outputs. The self-sustained nature of CPGs is believed to reduce 
the computation burden of the brain. As illustrated in Fig. 1, the central controller, similar to the brain of 
an animal, is expected to stabilize the vehicle dynamics by commanding a reduced number of variables such 
as the frequency {(jo{t)) and phase difference (A^j) of the oscillators instead of directly controlling multiple 
joints. The existence of CPGs has been conflrmed by biologists. ^^^^ Experiments with limbed vertebrates 
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have shown that individual hmbs can produce rhythmic movements endogenously.^^'"^^ Such empirical data 
have been interpreted as evidence that each limb has its own CPGs that can behave in a self-sustained way. 
However, sensory feedback is also known to play a crucial role in altering motor patterns^^'^^ to cope with 
environmental perturbations. Incorporation of sensory feedback into the CPG model has been presented in 
companion papers^^'^^ for a turtle robot. 

The most popular animal model for CPGs has been the lamprey, a primitive eel-like fish.^'' While the 
robotics community eagerly embraced the concept of CPG models for swimming or walking robots, ^^'^^^^ 
this work reports the first CPG-based control for flapping flight. The use of nonlinear oscillators for insect 
flapping flight has also been suggested by some biologist s.^^'^^ Clearly, flapping flight is technically more 
challenging to mimic than swimming and walking, due to its uncompromising aerodynamic characteristics. 



II.A. Mathematical Model of CPGs by Hopf Oscillators 

Our neurobiologically inspired approach centers on deriving an effective mathematical model of CPGs based 
on coupled nonlinear limit cycle dynamics. Once neurons form reciprocally inhibiting relations, they oscillate 
and spike periodically. Some prior work uses a discrete nonlinear equation that describes spiking and spiking- 
bursting of a neuron model. In the present paper, an abstract mathematical model of complicated neuron 
models is obtained by coupled nonlinear limit cycles that essentially exhibit the rhythmic behaviors of coupled 
neuronal networks. In the field of nonlinear dynamics, a limit cycle is defined as an isolated closed trajectory 
that exhibits self-sustained oscillation.^^' If stable, small perturbations (initial conditions) will be forgotten 
and the trajectories will converge to the limit cycle. This superior robustness makes a limit cycle an ideal 
simplified dynamic model of CPGs. 

In the present paper, we use the following limit-cycle model called the Hopf oscillator, named after the 
supercritical Hopf bifurcation model with a = 1: 




-A 



-A 



-co{t) 



Equivalently, x = f (x; ct) + u(t), with x = (u 




(1) 



where the limiting cycle of the Hopf oscillator is a circle with a radius p > 0, whose center is at (a, 0)-^. 
Also, A > denotes the convergence rate to the circular limit cycle and u(t) is an external or coupling input. 

For a single Hopf oscillator with u(t) = 0, a Lyapunov function V = can be used to prove 

global asymptotic stability to the circular limit cycle. Also, the bifurcation parameter a can switch from 

1 to — 1 such that This would change the stable limit cycle dynamics to the dynamics 

with a globally stable equilibrium point at the bias "a" (see [51]). Such a change can be used to turn the 
flapping oscillatory motion to the gliding mode, as shall be seen in Section IV. We assume (7 = 1 unless 
noted otherwise. For coupled Hopf oscillators, the stability proof is much more involved and discussed in 
Section II.C. 

The possibly time- varying parameter u;{t) > determines the oscillation frequency of the limit cycle. 
A time- varying a{t) sets the bias to the limit cycle such that it converges to u{t) = pcos {cot -\- S) -\- a and 
v{t) = psin {ujt + (5) on a circle. This bias "a" does not change the results of the stability proof. The output 
variable to generate the desired oscillatory motion of each joint is the first state u from the Hopf oscillator 
model in Eq. (1). 

There are two specific reasons why we prefer the Hopf oscillator to construct engineered CPG arrays 
(e.g., see the salamander robot"^^'^^ and the turtle robot"^^'^^). One essential property of the Hopf oscillator 
for our synchronization stability proof is that its limit cycle is a symmetric circle as opposed to Van der 
Pqj35 Rayleigh oscillators.^-^ Hence, we can verify the following: 



f (R(A)x; p; a) = R(A)f (x; p; a) R(A) 



cos A 
sin A 



— sin A 
cos A 



(2) 



where R(A) G SO{2) is a 2D rotational transformation such that R(-A) = R"H^) = ^^{^)- Additionally, 
this symmetry prevents the classical problem of instability due to switching between two stable systems, as 



4 of 24 



American Institute of Aeronautics and Astronautics 



one of the systems would need to be asymmetric.^ The stabihty of coupled Hopf oscillators has also been 
investigated in [45,46,55]. 

Second, our stability proof of phase synchronization exploits a scaling factor of a symmetric oscillator 

f (^x; p; cr) = gf (x; p/g; a) . (3) 
II.B. Rationales for CPG-Based Flapping Flight Control 

We articulate the three key advantages of CPG-based control for control of a flapping flying vehicle with 
complex wing kinematics: (1) reduced dimensionality, (2) adaptive pattern modulation, and (3) synchro- 
nization. 

II.B.l. Reduced Dimensionality and Bandwidth Requirement 

The CPGs in animal spinal cords are known to relieve the computation burden of locomotion in the brain. ^"^'^^ 
Similarly, one significant advantage of CPG-based control over conventional control approaches is that CPG- 
based control reduces the dimensionality and bandwidth of signals required from the main controller to its 
actuators. As shown in Fig. 1, the main outer-loop flight controller needs to command only the reduced 
number of CPG parameters (e.g., frequency, phase lag, and amplitude) and much less frequently, instead 
of directly commanding time-speciflc reference signals for all the degrees of freedom in the wings and the 
body. Hence, the use of engineered CPGs can be regarded as implementing motion primitives'^ for systems 
with many interacting degrees of freedom. For example, a learning-based controller'^ using CPGs needs to 
adapt only the reduced dimensional CPG parameters. Such a model reduction approach for flight control 
has not been exploited in the literature. The reduced dimensionality of the CPG-based approach makes 
learning-based adaptive flight control more practical. 

II.B. 2. Adaptive Pattern Modulation 

Birds and bats modulate the CPG parameters (frequency, phase difference, bias, and amplitude) for the 
flapping, twisting, lead-lag, cambering, and flexing of the wings during their flight, as a function of flight 
speed^''^''^ and flight modes (e.g., turning, gliding, soaring, cruising, hovering, and perching). High-speed 
fllm analyses^' reveal that the flapping angle and frequency are largest at zero forward speed or in hovering 
flight, and decrease with increasing flight speed V (e.g., oc some bats'^). Such time- varying 

CPG parameters, shown in Fig. 1, will change the shape, size, and flexing of the wings, which constitute 
the morphological flight parameters.' Prior studies in flapping flight, although true in steady flight, assume 
that there is a constant or very narrow range of optimal frequency or amplitude.^'^'^"-^^' -^^'^'''^^ Agile 
vehicles with multiple flight modes may require a large envelope or discontinuous parameter changes. Typical 
sinusoidal signals can be modulated as well for continuous and discontinuous changes in frequency by using 
an integral of frequency, but discontinuous changes in amplitude or bias would require a low-pass fllter, 
thereby adding phase delays. The Hopf oscillator guarantees continuous transitions for any time- varying 
input of these parameters, without a need for any low-pass flltering. This same behavior extends to the 
ability to handle any initial conditions and the ability to reject disturbances. 

II.B. 3. Synchronization of Symmetric and Symmetry- Breaking Oscillation 

Bats exhibit complex wing flapping motions generated by their multijointed and compliant wings. One aim 
of the neurobiological approach to engineered flapping flight is to produce the analytical model of a wing 
beat oscillator that matches empirical data.-^^' ■^''''^'^^ While the benefits of nonlinear limit cycles for CPG 
models are articulated above, deriving an effective CPG model for engineered fiapping fiight has been largely 
an open problem (e.g., limit cycle dynamics, network topology, and how to integrate input and feedback 
signals). The key research issues include how to ensure the amplitude or phase synchronization of multiple 
coupled CPG oscillators and how to opportunistically break the symmetry of the oscillators for performing 
agile maneuvering of a small fiapping fiight vehicle. Unique to this formulation of coupled oscillators is the 
ability to set and control the phase differences between multiple oscillators (see Aij in Fig. 1). In contrast, 
the use of an independently controlled sinusoidal function for each wing joint would not permit a control 
law that uses a phase difference A^j. Furthermore, stability proofs of coupled sine functions are largely open 
problems and more difficult to handle than the standard form of limit cycle oscillators such as x = f (x)+u(t). 
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Observations of birds have found these phase differences to be key in performing maneuvers, and we wih 
later see that they can be useful for vehicle stability. In the next section, we present how to construct stable 
coupled Hopf oscillators for the purpose of adaptive pattern modulation. 



II. C. Global Exponential Synchronization of CPG Oscillators 

Synchronization means an exact match of the scaled amplitude with a desired phase difference in this paper. 
Hence, phase synchronization permits different actuators to oscillate at the same frequency but with a 
prescribed phase lead or lag. However, a sinusoidal function is not adequate to entail the complex coupling 
and synchronization between various joints and limbs. Hence, the use of coupled nonlinear oscillators in this 
paper provides a feasible solution to construct complex synchronized motions of multiple wing joints. In 
essence, each CPG dynamic model in Eq. (1) is responsible for generating the limiting oscillatory behavior of 
a corresponding joint, and the diffusive coupling among CPGs reinforces phase synchronization. For example, 
the flapping angle has roughly a 90-degree phase difference with the pitching joint to maintain the positive 
angle of attack (see the actual data from birds in [3]). The oscillators are connected through diffusive 
couplings, and the i-th Hopf oscillator can be rewritten with a diffusive coupling with the phase-rotated 
neighbor. 

Pi- 



±i = f(xi;pi) - k{t) ^ 
jeAfi 



X, - f R(A,,)x, 



(4) 



where the Hopf oscillator dynamics f (x^; pi) with a = 1 is defined in Eq. (1), A/^ denotes the set that contains 
only the local neighbors of the i-th Hopf oscillator, and rrii is the number of the neighbors. The 2x2 matrix 
R(Aij) is a 2-D rotational transformation of the phase difference Aij between the i-th and j-th oscillators. 
The positive (or negative) Aij indicates how much phase the i-th member leads (or lags) from the j-th 



member and A.- 



The positive scalar k{t) denotes the coupling gain and can be time- varying for 



different flight modes. 

We construct as many degrees of freedom as needed to more accurately model the joints of the wings, 
but let us focus on the key three flapping motions deflned in Fig. 2, namely flapping angle 0^, wing pitch 
(twisting) angle 0^^ and wing lead-lag angle i/jw Additionally, we assume that there is a second flapping 
joint (I)yj2 in the wing that can reduce the drag in the upstroke by folding the wings toward the body. Then, 
we can construct the whole state vector of the coupled oscillator such as 



{x} = 













/ {(t)yj^ -ai,vi)^\ 


X2 










{OwR - 02,V2Y 


X3 




{U3 






(V^^^ - as, vs)^ 


X4 




{U4 


- o^.vaY 




{(t>w2R - 04.,V4.Y 


X5 












X6 




{ue 


- Oq.VqY 




{Owl - Oq,Vq)^ 


X7 






- 07, Vj)^ 



















(5) 



Note that x^ here might represent the shifted Hopf oscillator vector such that x^ = (ui — ai^vi)^ as seen 
in Eq. (1), where ai{t) is the center of oscillation. For example, if we need a 10-degree offset for the main 
flapping stroke angle (jy^^ then we can set ai = as = 10 deg. so that the flapping stroke angle oscillate around 
10 degrees. 

In order to analyze the stability of a network, we need to construct fully coupled dynamics of the 
augmented state vector {x} from Eq. (4) 



{x} = [f({x};p)]-A:(t)G{x} 



(6) 



where [f({x};/))] = [f (xi; pi); f (x2; P2); • • • ; f (x^; p^)]. Also, p = {pi,-- - , pnY ^^"^ implying that the radius 
of each oscillator can be different. The 2n x 2n matrix G is a Laplacian matrix with phase shifts R(A^j) 
constructed from Eq. (4). 

The coupling topology and phase shift between each oscillators are reflected in the G matrix. Such 
phase shifts along with the bifurcation parameter a can be used to deflne different flight modes, similar 
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(a) Configuration A (b) Symmetric Configuration A (c) Configuration B 



Figure 3. Graph configurations of the coupled Hopf oscillators on balanced graphs. Many other configurations are 
permitted in this paper and the unidirectional couplings can be replaced by the bi-directional couplings. The numbers 
next to the arrows indicate the phase shift Aij from the i-th member to the j-th member while Figure b shows 
the nominal values of the phase shift from the symmetric wing configuration such that A21 = Ass =90 deg. and 
A31 = A75 = —90 deg. Such phase shifts define flight modes (wing movement gaits). 

to walking gaits. Numerous configurations are possible as long as they are on balanced graphs^ and we 
can choose either a bidirectional or a uni-directional coupling between the oscillators. Some configurations 
considered in this paper are shown in Fig. 3. The numbers next to the arrows indicate the phase shift A^j, 
hence Aij > indicates how much phase the i-th member leads. Since the graphs in Fig. 3 are on balanced 
graphs, the number of input ports equal the number of output ports. Further, all the phase shifts (A^j) 
along one cycle should add up to a modulo of 27r. Figure 3b shows the nominal values of the phase shift from 
the symmetric wing configuration such that A21 = Ags = 90 deg. and A31 = A75 = —90 deg. The empirical 
data suggest that the pitching angle {O^j) has approximately a 90-degree phase lag with the flapping angle 
(^^y), which agrees with the aerodynamically optimal value. ^'^^ For hovering flight, Dickison,^^ using his 
Robofly testbed and numerical simulations, found that increasing the phase difference value A21 to 90 deg -\-S 
further contributed to enhancing the lift generation, which is explained by the wake capture and rotational 
circulation lift mechanism. Hence, the ability to control A21 allows us to investigate the optimal value of the 
phase difference. In addition, the nominal value of A31 = —90 deg, the phase difference between the flapping 
stroke angle and lead-lag angle will results an elliptical orbit of the wing. On the other hand, by having two 
difference phase differences for the left and right wings, we can investigate how symmetric-breaking wing 
rotations contribute the agile turning of flapping flight. Furthermore, by having an independent control of 
the phase difference A31 and A75, we can investigate another symmetry-breaking impact of the differential 
delay in the lead-lag motion. Such differential phases can be used to stabilize the flapping flying dynamics. 
The G matrix in Eq. (6) for Fig. 3a can be found as G = 



2I2 








^R(A3i) 















I2 























-gR(A3i-A2i) 


I2 























P3 ^ 


I2 














pi ^ 











2I2 
























I2 























-^R(Ar5-A65) 


I2 























P7 ^ 


I2 














(7) 



where in general, the radii (the amplitude of the oscillation from the bias a^) are symmetric such that pi = p2^ 
p2 = pA = p7i and p5 = although the difference of the maximum amplitude of each oscillation can 
occasionally be used to generate side forces or turning (rolling or yawing) moments. 

The proof of phase synchronization boils down to finding the condition of k by which the flow- invariant 
synchronized state, constructed from G{x} = 0, is globally stable. In fact, by using contraction theory, ^^'^^ 
we can prove global exponential synchronization of the coupled Hopf oscillators. We first introduce the main 
theorem of contraction theory. 

Theorem 1 For the system x = f(x, t)^ if there exists a uniformly positive definite metric, M(x, t) = 
0(x, t)^0(x, t)^ where is some smooth coordinate transformation of the virtual displacement, 87^ = 0(5x^ 
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such that the associated generalized Jacohian, F is uniformly negative definite, i.e., 3£ > such that 



0(x,i) + 0(x,i)— ) 0(x,i) '<-a. 



(8) 



then all system trajectories converge globally to a single trajectory exponentially fast regardless of the initial 
conditions, with a global exponential convergence rate of the largest eigenvalues of the symmetric part ofF. 

Such a system is said to be contracting. 

Proof 1 The proof is given in [62] by computing ^Sz^Sz = 2Sz^FSz. □ 

The synchronized flow-invariant subspace for the conflguration in Fig 3a is deflned by G{x} = such 



that 



X({x}) 



Pi 

P2 



Pi 



Pi- 



R(A12)X2 = ^R(A13)X3 = ^R(A13)X4 



P3 



P4 



X5 = ^R(A56)X6 = ^R(A57)xr = ^R(A57)X8 



pi 



pi- 



P6 



P7 



(9) 



where we used Aij = —Aji. 

The flow invariant subspace M. in Eq. (9) can be re-written with respect to the flrst state vector xi = zi 
such that 

X({x})^zi=Z2 = --- = z,, {z} = T(A,,-,pO{x} (10) 

^R(Ai3)x3, and so on. For example, 



where {z} = (zi, Z2, • • • , z^)^ and zi = xi, Z2 = ^R(Ai2)x2, Z3 



the T matrix for the conflguration in Fig. 3a is given as 
T{Aij,pi) = diag 



I2, ^R(Ai2), ^R(Ai3), ^R(Ai3), ^l2, ^R(A56), ^R(A57), ^R(A57) 
P2 P3 P4 Pb P6 P7 P8 



(11) 



Then, we present the main theorem of this section. 

Theorem 2 // the following condition is met uniformly for \/t > 0, any initial condition {x} of the coupled 
Hopf oscillators in Eq. (4) and Eq. (6) on a balanced graph converges to the flow-invariant synchronized 
state Ai exponentially fast. 

mXmin (V^(L + L^)V/2) > A (12) 

where A is the convergence rate of the Hopf oscillator in Eq. (1), Xmin (V-^(L + L-^)V/2) denotes the 
minimum eigenvalue, and L is the Laplacian matrix constructed from the balanced graph such that G = 
T~^LT withT defined from Eq. (10). In addition, the real orthonormal2nx2{n—l) matrixY is constructed 
from the orthonormal eigenvectors of (L + L'^)/2 other than the ones vector 1 = (I2;l2; • • • ;l2) such that 
VV^ + llVn = l2„. 

Proof 2 The proof can be obtained by exploiting results from [55]. The proof here is simpler than [45, 4^] in 
the sense that we derive the Laplacian matrix and orthonormal flow-invariant matrix that are independent of 
the rotational angles. Consider the orthonormal space W , constructed from the orthornomal eigenvectors of 
the symmetric part ofL (see [6]). Then, the global exponential convergence to the flow-invariant synchronized 
state Ai is equivalent to 



V {z} 0, globally and exponentially 
By pre- multiplying Eq. (6) by T and using T{x} = {z} and G = T^-'^LT^ we can obtain 

{z} = T[f({x};p)]-fc(i)L{z} 

where the CPG network in the example in Fig. 3a is on a balanced graph such that 
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In other words, we transformed the G matrix to the conventional graph Laplacian matrix L. 

We recall the definition [f({x};p)] = [f (xi; pi); f (x2; ^2); • • • ;f(xn;Pn)]- 
Since T [f ({x}; p)] = T [f (T-i{z}; p)] , we can find 



T[f({x};p)] 



Pi 



R(-Ai,-)f(x,;pO 



^R(-Ai,)f(^R(Ai,)z,;pO 
Pi pi 



(16) 



[f ({z}; pil)] = [f (zi; pi); f (zs; pi); • • • ; f (z^; pi)] 



where we used f(R(A)x) = R(A)f(x) and f (^x; p) = ^f(x; p/g) from Eqs. (2) and (3). 

It is very important to notice that the radii of the final augmented Hopf oscillators in Eq. (16) are 
identical to pi such that p = pil. Such a transformation to an identical oscillator of the same radius pi is 
indispensable for the following proof of synchronization. Also, Eq. (16) can he obtained only by exploiting 
the circular symmetry of the Hopf limit cycle, i.e., Eqs. (2) and (3). 

By pre-multiplying Eq. (I4) with V-^ and substituting {z} = VV-^{z} + ll-^{z} into Eq. (16), we can obtain 



V^{z} = [f(VV^{z} + ll^/n{z};pil)] - A:(t)V^LVV^{z} 

where we used Lll^ = 0. 

We can construct the following virtual dynamics of y from the preceding equation 

y = [f(Vy + ll^/n{z};pil)] - k{t)\^L\y 



(17) 



(18) 



which has y = \^{z} as one particular solution. 

The other particular solution of Eq. (18) is y = since [f(ll^/n{z};pil)] = 0. If the virtual 
system y in Eq. (18) is contracting, then the two particular solutions tend to each other exponentially from 
any initial condition, thereby completing the proof. 

The virtual system Eq. (18) is contracting ifY^ [f] V - A:V^(L + L^)V/2 < due to Theorem 1. This 
condition is equivalent to kXmin (^"^(L + L-^)V/2) > since the maximum eigenvalue XmaxC^^ [f] V) < A. 
For the example network in Fig. 3a, this condition corresponds to k > A/0.198. This also holds for a time- 
varying k(t), if the condition holds uniformly for all time. 

The same proof works for an arbitrary CPG network on a balanced graph that has V^(L + L^)V/2 > 0. 
For undirected graphs (all the connections are bi-directional), L automatically becomes a balanced symmetric 
matrix. □ 



In conclusion, Theorem 2 can be used to find the proper couphng strength k to exponentiaUy and globaUy 
stabihze the coupled Hopf oscillators given in Eq. (4). 



II. D. Fast Inhibition of Oscillation by Hopf Bifurcation 

As stated earlier, we can rapidly inhibit the oscillatory motion of the coupled Hopf oscillators in Eq. (4) by 
exploiting the bifurcation property of the Hopf oscillator model. In other words, changing the <j = 1 in Eq. 
(1) to a = —1 would rapidly convert the limit cycle dynamics to exponentially stable dynamics converging 
to the origin such that u ^ a and v ^ This single bifurcation parameter (cr) can be used to switch the 
flapping flight mode to the gliding mode or soaring mode without dramatically changing the CPG oscillator 
network. Simulation results that alternate between two different flight modes are presented in Sec. IV. 

Theorem 3 For any positive gain k > 0, any initial condition {x} of the coupled Hopf model with a = —1 
given in Eq. (4) converges to the origin ({x.} 0) such that Ui ai and Vi ^ for a// i = 1, • • • , n. The 
oscillation frequency uj need not change to zero. 



Proof 3 It is straightforward to show that a = —1 will make the uncoupled Hopf oscillator in Eq. (1) 

exponentially stable dynamics for any {u, v) except the shifted origin (a, 0) since the symmetric part of the 
Jacobian F in Eq. (8) is now strictly negative definite regardless of any uo. Thus, any positive k will lead 
to exponentially synchronizing dynamics that tend exponentially to the shifted origin (a^ , 0) , and this can be 
shown similar to the proof of Theorem 2. □ 
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(a) 



(b) 



(c) 



Figure 4. A 8-DOF robotic bat (with 10 control variables) MAV' 



20 



We can also turn the limit cycle dynamics to the dynamics with a stable equilibrium by changing the 
coupling gains, as described as fast inhibition in [6]. However, the method using bifurcation is superior in 
the sense that we can keep the original coupling gains and Laplacian matrices for alternating flight modes. 
It should be noted that changing uo to zero would also result in no reciprocal flapping motion, however the 
converged steady-state value depends on the initial conditions, whereas a = — 1 would always make the 
system to converge to the desired value (a^,0) for each joint i. 

II. E. Perspectives on Sensory Feedback Connection 

The property of robustness, inherent in the CPG-based control, is particularly emphasized by the literature 
(see [63]). Stable locomotion can be achieved using the interaction between the CFG model, the physical 
model of the body, and the environment.^^ Most models^^'^^'^^ use an open-loop approach without sensor 
feedback, while some others^^'^^ incorporate sensor feedback to modulate the reference oscillator patterns. 
One drawback is that such open-loop approaches do not ensure the synchronization of the physical states 
in the presence of external disturbances. In other words, the mutual entrainment^^' between the CPG 
and the mechanical body is not guaranteed. Recently, two new closed-loop CPG control strategies that 
reinforce emerging rhythmic patterns of actual physical joints like foil-fin actuators has been proposed. ^^'^^ 
Such a reflex-based closed-loop CPG method has a potential for discovering practical ways of flapping wing 
coordination in the presence of external disturbances. For example, we can modify Eq. (4) to incorporate 
actual states of actuators: 



where and pj are the actual actuator joint angles measured directly from sensors and Vi/vj denotes the 
new amplitude scaling for actuators. Then, similar to Theorem 2, we can prove the stability of Eq. (19) 
by assuming the actuator states can be represented by p^ = p^(t)R(A^^iag)x^ where A^^iag indicates a time 
delay between p^ and x^. Readers are referred to [45,46] for details. 

In Section III, we present the wing kinematic model and the dynamic model of flapping flight dynamics 
that can be driven by the CPG network. 

III. Wing Kinematics, Aerodynamic Forces, and Vehicle Dynamics 

We derive a three dimensional model of the wing kinematics in this section. The wing kinematic model 
supports both flexible and rigid wings. Also, we present the 6-DOF dynamic equations of motion of flapping 
flight that can be used to validate the coupled wing control driven by CPG. We also show that the effective 
angle of attack of each blade element can be effectively controlled by the synchronized pitching (6>^) control 
as well as the wing joint angles such as (p^ and ipw. 

We present a realistic modeling that encompasses a tilted stroke angle, the lead-lag motion, and the 
relative body velocity, in addition to the stroke and pitch angles. In deriving these equations, the actual 
control degrees-of-freedom of the robotic bat MAV testbed shown in Fig. 4 are considered. Its half wingspan 
is 32 cm. The 8-DOF robotic bat has 10 independent control variables including separate flapping frequency 
and amplitude control for each wing as well as pitch and lead-lag angle servos (see [20] for details). 




(19) 



10 of 24 



American Institute of Aeronautics and Astronautics 



Figure 5a shows a side view of the flapping flying MAV with the body frame X5 = 7/5, ^5)-^ and the 
stroke-plane frame = {xg^Vs^ ^s)^ of the right wing. Similar to the robotic bat in Fig. 4, we assume that 
each wing has flapping {4>w), pitch (0^)^ and lead-lag (i/j^) control. We develop only the equations of the 
right wing since the similar expressions for the left wing can straightforwardly follow. The center of the 
stroke-plane frame is located at {dx^dy^dz), and it is tilted by the inclination angle ©s(0, which can be a 
function of time and the forward velocity. Without the lead-lag motion, the axes i/s and Zs define the stroke 
plane. Hence, the transformation between these coordinate axes can be given by 



Xfe = Tbs{Qs)^s + {dx, dy, dzY , where Tiys{Qs) 



cos Qs sin 6^ 

1 
- sin © s cos © c 



(20) 



where in this paper denotes the transformation from x^ to x^, whereas T^^ = T£ would correspond to 
the transformation from xj, to x^. 

For a hovering insect, the stroke plane is almost horizontal (i.e., ©^ = 90 deg in our coordinate definition 
in Fig 5a), resulting in forward and backward reciprocating motions. This is the assumption used for some 
prior work.'''^^'-^^'-^^ In contrast, the stroke angle of birds and bats varies as a function of flight speed; at a 
low speed, the angle is almost horizontal (©^ = 90 deg) and it approaches ©^ = deg as the flight speed 
increases. 

If there is no lead-lag motion, the additional transformation for a wing stroke angle would complete 
all the required transformation between the body frame and the wing frame. However, a nonzero lead-lag 
angle further complicates the wing kinematics. Choosing the rotational axes for flapping, lead-lag, and 
pitch depends on the actual hardware setup and actuators, and our choice is influenced by the robotic bat 
MAV shown in Fig. 4 (see [20]). In contrast with Azuma's derivation in [3] where the stroke angle ©s(t) 
is dependent on the (j)w{t) and the lead-lag angle ipwit)^ our ©s(t) is an independent control variable. Our 
decision is based on the observation that ©s(t) can be an important control variable for eflicient engineered 
flapping flight. Further, this kind of actuator mechanism is easier to implement and control. As shown in 
Fig. 5b, the lead-lag angle is deflned by the rotation about the Zs axis- the z-axis in the stroke plane frame. 
In contrast with the flxed angle rotation in [3] , then we rotate about the new x-axis to obtain the wing frame 
x^. For both wings, the positive direction of is the forward direction, while the positive stroke angle 
(pw indicates an upstroke motion. This sign convention does not agree with the original positive direction of 
rotation for the right wing, so extra care should be taken to determine the correct angular transformation 
matrices. 

For the right wing, the transformation between the stroke plane frame (x^) and the wing frame (x^) can 
be written as 



COS sin 

- sin 7/; COS 
1 



1 

cos^ 
— sin ^ 





sin^ 
cos d 



(21) 



In order to compute the local lift and drag of a blade element with width dr and wingspan coordinate 
r e [0,i?], we need to transform the velocities in body coordinates to the incident velocities in the rotated 
wing frame. For example, consider the vehicle body speed Vb with the body angle of attack ax and the 
side-slip angle ay. Note that ay is commonly denoted by (3 in the aerospace community, but in this paper 
(3 denotes the direction of the relative wind of a blade element. Then, the free-stream velocity in the body 
frame can be written as 



= {Vb cos ay cos ax, Vb sin ay, cos ay sina^^)^ + + 



(22) 



where and v^; denote the induced velocity vector and the wind velocity vector respectively. In other 
words, in the absence of and v^;, the vector V5 equals the velocity of the vehicle in the body frame. Let 
us assume that ax and ay include the effects of the induced velocity and ve is small. 

Then, the free- velocity vector V5 in the body frame can be transformed to the wing frame. In addition, we 
can also compute the additional velocity on the wing frame induced from the body angular rate ft^, = (p, r)^ 
and the offset distance d = (d^, dy, dz)^ of the stroke plane frame (see Fig. 5a). By adding these two terms, 
we can obtain 

= T^,((/)^, V^^)T,5(©,)(V5 + X d) (23) 
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(a) Transformation from the vehicle body frame to the stroke plane frame of the right wing. 




(b) Transformation from the stroke plane frame to the wing frame . 

Figure 5. Schematic of the 3D wing motions 
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In order to compute the rotational velocity on the wing frame produced by the flapping 0^ and lead-lag 
ipw motions, as well as a relatively slower stroke angle change ©s(t), it is more convenient to construct the 
angular rate vector in the stroke plane frame as follows 



-cos llJw(pw 

^tot = T^sb{Qs)^b + I sixiip^ij)^ + 6s 



(24) 



Then, we can compute the induced rotational velocity from the wing motions of the blade element dr 




(25) 



where ^^^(r), ywi^i^) and Zuj{r) are the deformation of the blade element due to aeroelastic deformation or 
active cambering control that can be found in bat flight. Hence, the derivations in this section can be used 
for flexible wing models, although the Cl{o^) and Cd{<^) functions should be corrected for such cambered 
wing shapes. 

By adding in Eq. (23) and V^^^ in Eq. (25), we can obtain the total velocity of the wind at the 
blade element, located at r on the wing span axis. 



V — 




= T^s{(t>w.i^w)Tsb{Qs)(Vb^^bXd)^{T^s{(t>w.^w)^tot) 




ijwir) 
(26) 

A similar expressions can be obtained for the left wing. Also, the preceding derivations can straightforwardly 
be extended to account for the second joint (elbow) if each wing has one. 

Now, we can obtain the local effective angle of attack ayj of the blade element to determine aerodynamic 
forces and torque. Let us assume that the deformation of a rigid wing is negligible and there is no active 
cambering control. Also, the contribution from the body angular rate il^, is small. Then, Eq. (26) reduces 
to 




= T^s{(t>w,'^w)T^sb{QsWb + 



^ / -cosV^^^^ 

ws 







f1 




X 


r 


1) 







(27) 



Then, we can obtain the local incident angle (3^, (measured clockwise), the angle of attack a^, and the speed 
of the wind Vr on the blade element on the right wing as follows 



(3w{r,t) = tan 



aw{r,t) = 0w{t) - (3w{r,t) 



(28) 



where we neglected the flow along the wing span V^y and the wing rotation 0^ {t) controller can be properly 
designed to yield a positive angle of attack for both upstroke and downstroke motions. 
If the MAV were flying with zero flight path angle and 6^ = '0^ = 0, we can obtain 



, -1 rcpw 

= tan — 



- tan 



and kr = — :- 



(29) 



where the reduced frequency kr compares the velocity by the wing flapping motion with the forward speed, 
thereby indicating the degree of unsteady aerodynamics (if kr ^ 1, unsteady effects dominate^' ^^). 

We can see that the sign of (3^ is consistent with the positive direction of the flapping (stroke) angle (p^ 
since the downstroke < leads to the negative flow angle P^j < 0- Also, becomes a function of 0^ 
indicating why the pitch angle 0^ should have a 90-deg phase difference with (j)^ (A21 =90 deg) in order to 
maximize the local angle of attack. Figure 6 shows a result of the wind-tunnel test of the robotic bat shown 
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in Fig. 4. The robotic bat is mounted horizontally in a wind tunnel [20]. In Fig. 6, the pitch motion {0^) 
was activated in 30 sec. with about A21 =90deg. This increased the lift (F^) by more than a factor of two 
while the wind speed and the CPG frequency u were held constant. This result shows that the synchronized 
pitch control is indispensable. We can find the most efficient phase difference (A21), which can be compared 
with Dickinson's experimentation with robotffy (advanced rotation with A21 = 90 -\- 5 deg). 




20 25 
Time [s] 



Figure 6. Impact of the synchronized pitch 9u, oscillation. 



Once we obtain the local effective angle of attack , we can proceed to obtain the aerodynamic forces of 
the blade element by evaluating the lift and drag coefficients, Cl{(^w) and Cd{o^w)' Flapping ffight, typically 
within a low Reynolds number regime (Re < 10^), is governed by unsteady aerodynamics characterized by 
large-scale vortex structures. It is understood that the main lift enhancement mechanism of ffapping ffight is 
governed by (1) the leading edge vortex (LEV) that leads to delayed stall at a very high angle of attack, (2) 
the rotational circulation lift, and (3) wake capture that generate aerodynamic forces during ffapping angle 
reversals. In particular, Dickinson's series of papers,^' by cross- validating the numerical computation and 
experimentation using the Roboffy, shows that a quasi-steady aerodynamic model predicts the aerodynamic 
coefficients reasonably well. While computational Fluid Dynamics (CFD) is much too computationally 
burdensome to justify its use in control design, this quasi-steady approximation method can be veriffed and 
improved by the experimental set-up such as the robotic bat described in [20]. 

The seminal paper by Dickinson^^ used a hovering pair of wings without a forward speed as follows 

Cl{(^u,) = 0.225 + 1.58 sin(2.13a^ - 7.2 deg) 

Coiaw) = 1.92 - 1.55cos(2.04a^ - 9.82 deg) ^ ^ 

It should be noted that Dickinson's robotffy's setup used a horizontal stroke plane, as typically seen in insect 
ffight, whereas we assume a 90-deg stroke plane angle. The angle for a general ffapping wing is time- 
varying, as described in this section. A recent paper^ considers a nonzero-forward speed and the coefficients 
of Eq. (30) can be modiffed to become functions of the reduced frequency (/c^). 

From the quasi-steady approximation, we can compute the lift and drag forces acting on the blade element 
with width dr as follows. 

dL = ^pCl (aUr.t)) c{r)V,\r,t)dr, dD = ^pCo (aUr^t)) c{r)V,\r,t)dr (31) 

In addition, Ellington^ derived the wing circulation F^ = 7rdc^(3/4 — xq) based on the Kutta-Joukowski 
condition. This quasi-steady approximation for the rotational lift can be written as 

dLrot = (^27r(^ - fo)^ c^{r)Vr{r,t)au:dr (32) 

where xq is the location of the pitch axis along the mean chord length. Also, can be computed from Eq. 
(28) and often approximated reasonably well by the angular rate of the wing pitch motion 0^. 

The total x and z directional forces of a single wing (either right or left) in the body frame are obtained 

as 

nR pR 

Fwz = / dD sin - {dL + dLrot) cosPuj, F^jx = / -{dL + dLrot) sinPuj - dDcosP^ (33) 
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where the positive direction of is downward as shown in Fig. 5a. 

The F^x and F^^ forces on the wing frame given in Eq. (33) can be transformed into the forces in the 
vehicle body frame: 



'Fx' 



right 



right 







(34) 



right 



where we added the subscript right to indicate the right wing. A similar expression can be obtained for the 
left wing (Fieft). Each wing has different wing angular parameters such as (j)^^ V^^, and 0^^ although the 
stroke plane angle O5 is the same for both wings. 

In order to compute the rotational moments generated by the aerodynamic forces, we first calculate the 
position of the wing blade element with respect to the body frame 



p(r) Tbs{Qs)Tsw{(l)w^^w) 



Then, we can compute the aerodynamic moments with respect to the e.g. 




(35) 



dMy 

(dM^o\ 

dMyo 
XdM^I 



p(r) X Tk,{Q,)T,^{(t)y„%l)^) 



-{dL + dLrot) — dD cos \ 



dD sin/3u, — {dL + dLrot) cos / 
( rCio 

rCno 



dMyo 



f 



dMy, 



f 

Jr=0 



\ 



dM, 



(36) 



(37) 



(38) 



where dMxo, dMxo, and dMxo denote the constant aerodynamic moments that include the moment at 
the mean aerodynamic center, computed by the moment coefficients C/o, Cmo^ Cma,w^ and Cno- The 
transformation matrix Tq^{0^) rotates the wing frame about the axis by the wing pitch rotation angle 



By combining all the forces and moments from the right wing and the left wing, we can derive 6-DOF 
equations of motion for the flapping flying MAV in the body frame, whose orientation with respect to the 
inertial frame is described by the Euler angles. We assume the mass and the moment of inertia of the wing 
compared to the body weight are negligible so that the e.g. remains fixed. The translational motion of the 
e.g. of the vehicle driven by the aerodynamic force terms in Eq. (34) can be expressed as 



mVfe + X Vfe = Tfee(0fe, Ob, ipb)Fg + F^ght + Fieft + A 



(39) 



where all the symbols are defined in the nomenclature section. Each wing has different wing angular param- 
eters such as V^w, and 0^, although the stroke plane angle 65 is the same for each wing. 

The equations of rotational motion are driven by the aerodynamic moments M^ight and Mieft of each 
wing from Eq. (36) 

IbUb + 1^6 X (Ifel^fe) = Mright + Mieft + B (40) 

The relationship between the body angular rate fib = {p, g, r)^ and the Euler angle vector = (0^, ^5, ipb)'^ 
can be determined by^^ 



qb = z;(q6)ri6 



(41) 



where any other orientation representations such as quaternions can be used in lieu of the Euler angles in 
the preceding equations. Also, any disturbance force and torque can be added to the equations. 
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IV. CPG-based Flapping Flight Control and Simulation Results 



The aim of this section is to show that CPG-based flight control can stabilize and control flapping flight 
dynamics given in Sec. Ill by commanding a reduced set of CPG parameters that generate the phase- 
synchronized or symmetry-breaking oscillatory motions of two main wings. In particular, we show that the 
dynamics can be controlled without a tail or aerodynamic control surfaces such as ailerons, elevators, rudders, 
and directional control of tail wings. Another important contribution is to demonstrate stable transitions 
between gliding and flapping flight modes by using the Hopf bifurcation. 

IV.A. CPG-based Flapping Flight Control 

The example presented in this paper has three different flight modes: gliding, flapping, and turning. We 
show that only three control laws for uj{t), the symmetric phase difference (A32 = Ayg) between Oyj and ip^, 
and symmetric breaking of the lead-lag max amplitudes {ps 7^ pr) can control both longitudinal and lateral 
dynamics of tailless flapping flight with six independent wing joint angles. 

IV.A.l. Flapping Flight Control by Flapping Frequency: 

We propose a novel control law unique to our CPG set-up which reduces control dimensionality to only 
three parameters. The first parameter is the oscillation frequency u;{t) of the coupled Hopf oscillators in Eq. 
(4). The flapping frequency, correlates with increased lift and thrust. Those in turn correlate with the 
velocity of the body. For example, we can consider the following control law 

U;{t) = L0dt = K^ I (y^,desired " ^x,actual) dt (42) 

Jo Jo 

and use the following corollary. 

Corollary 1 From the dynamic equation of the Hopf oscillator in Eqs. (1) and (4), any time-varying ijj{t) 
does not affect the synchronization stability proof for Theorem 2. 

Proof 4 Since the symmetric part off cancels the uj{t) term and uo{t) does not change the maximum eigen- 
value of\^ [f] V. The rest of the proof follows Theorem 2 □. 

IV.A. 2. Flapping Flight Control by Phase Differences and Symmetry Breaking: 

The phase differences (A^j) of the coupled Hopf oscillators in Eq. (4) can be effective control mechanisms 
for flapping flight dynamics, and this controllability is another reason that justifies our coupled nonlinear 
oscillator framework rather than individual control of each angle. However, our synchronization stability 
proof in Theorem 2 assumes constant or relatively slowly varying phase differences. Our previous work^^ 
showed that the synchronization error terms from time- varying ^ij{t) are bounded and can be made arbi- 
trarily small with a sufficiently large coupling gain k due to robust contraction analysis. In the present 
paper, we show that the synchronization errors due to time- varying A^j(t) and pi{t) can tend exponentially 
to zero, thereby further strengthening our claim with the coupled Hopf oscillators. 

Corollary 2 Let us assume that the synchronization condition of k in Theorem 2 holds. For time-varying 
phase differences A^j(t) or time-varying radii Pi{t), pjit), the synchronization errors of the rotated Hopf states 
{z} globally and exponentially converge to zero, i.e., V-^{z} 0, if we add the additional time-varying term 
_rp-irp|^| coupled Hopf oscillator dynamics in Eq. (6) such that 

{x} = [f ({x}; p)] - fc(t)G{x} - T-it{x} (43) 

where the block diagonal transformation matrix T{Aij^pi) is defined such that {z} = T{x} as in Eq. (11). 
Also, T = ■^T{Aij, Pi) is a function of Aij and ^Pi/pj for each time-varying Aij, pi, and pj. 

As shall be seen later, we do not have to directly compute the time-derivatives of Aij{t) and pi{t)/pj{t) if 
we select the control laws carefully. 
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Proof 5 Since T{x} = {z}^ we can verify {z} = T{x} + T{x}. Hence, pre- multiplying Eq. (43) by T 
results in the same dynamics of {z} in Eq. (I4) even with time-varying Aij{t), pi{t), and Pj{t). The rest of 
the proof is identical to Theorem 2. □ 

If we select graph configuration A in Fig. 3, we have four symmetric phase differences available. We 
use two degrees of freedom as follows. We do not break symmetry of the right and left wings in the phase 
differences between the pitch and the lead-lag angles by setting 

A76 = A75 - A65 = A32 = A31 - A21 = -KAs^Ob + Ao (44) 

where Ob denotes the pitch Euler angle, i^Aas the control gain, and Aq is a nominal value that determines 
the steady-state pitch angle of the vehicle body (^^,). Not only is this effective in stabilizing the longitudinal 
motion, the nominal value can be used to select ascent or decent angle. 

In order to achieve perfect synchronization, we modify the coupled Hopf oscillators as discussed in 
Corollary 2. It is important to note that we do not have to differentiate A76 = A32 to obtain R(A32) in T. 
This follows from Eq. (44) 

Ar6{t) = A32{t) = KAs2^b = -Kas2 (<?cos^5 - r sin ^6) (45) 

where we used the relationship between the body angular rate and the rate of the Euler angles given in Eq. 
(41). 

In order to initiate a turning maneuver, consider breaking the symmetry of the right and left wings. One 
method of symmetry breaking is the maximum amplitude of the lead-lag angle (^3 and for the example 
in Fig. 3). We set this proportional to yaw rate to provide yaw rate damping. Due to the coupled nature of 
the motions, this causes both the yaw rate and the roll angle to go to zero. If we want to control the bank 
angle directly by setting desired ^5, we can use 

Psit) = P3, nominal " Kr{(t)h — 06,desired), Plif) = P7,nominal + ^r(06 " </>6,desired) (46) 

where ^6, desired cau be set to zero to end turning. The coupled Hopf oscillators again can be modified 
as in Corollary 2 to achieve zero synchronization errors with the time-varying radii. Also, we need not 
differentiate p^it) and pj{t) directly since <j)i) = p -\- g sin 0^ tan ^5 + r cos 0^ tan ^5. We used the method of 
symmetry breaking for our simulation. 

An alternative method of symmetry breaking is to add a difference value S to the nominal value (90 deg) 
of the phase differences of the flapping and pitch angles such that Aes = 90deg + ^ and A21 = 90deg — ^. The 
phase difference between flapping and pitch is vital for lift and thrust generation. Therefore, this difference 
between the right and left wings causes roll and proverse yaw. This method of symmetry breaking was 
observed in [61]. 

IV. A. 3. Gliding Mode Control 

We produce gliding flight with no reciprocal flapping motion by setting the bifurcation parameter cr = — 1 
in Eq. (1). As discussed in Sec. II. D, simply setting uj{t) =0 without a = —1 will not ensure convergence to 
the controlled bias values (a^). This provides us with simple control of our wing by exploiting the bifurcation 
of Hopf oscillators, causing them to exponentially converge to a single non-oscillatory value corresponding 
to the bias (a^). We can then control the lead lag motion (ipw), flapping angle {(pw), ^ind wing pitch angle 
(0) by their bias parameters. In gliding flight, synchronization is less critical so we can turn off the coupling 
by A: = 0. A negative (positive) flapping angle or negative (positive) lead-lag angle can provide a pitch-down 
(pitch- up) moment due to drag or lift, respectively. We have therefore reduced control dimensionality to 
three actively controlled parameters: wing pitch, wing flapping angle, and lead-lag angle. In fact, depending 
on the physical characteristics of the specific vehicle, controlling only one of wing flapping angle or lead-lag 
angle could be sufficient for longitudinal gliding stability. For the network in Fig. 3a, we can use the following 
bias for the lead-lag angle of each wing: 

asit) = o.7{t) = -kpOb - kdQ - ki Obdt + tpMas (47) 

Jo 

where /cp, kd^ and ki are positive gains and ^lJb^as is a constant. The simulation uses PID control of flapping 
angle and lead-lag bias. The angle which implies a dihedral angle in gliding flight, turns out to be the 
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most important, which is in keeping with a previous study of the importance of drag-based stabihty,^^ and 
the conclusion that birds ghde more hke taihess vehicles than conventional tailed aircraft. Only integral 
control is used for pitch bias, in order to obtain a constant wing angle of attack. 

IV. A. 4' Two Alternating Flight Modes for Altitude Control 

Inspired by altitude stabilization of animal flight, we propose a switching logic between flapping mode and 
gliding mode. Requirements for switching use the current bifurcation parameter a G (—1,1) to determine 
what mode we are in, as well as altitude and velocity information to determine whether to switch mode. 
Recalling that the z-direction is positive downward, we set the test for gliding mode as 

if (cr = 1 & Zb < -/lmax,flap & Vbx > V^x,max) Or ((7 = -1 & Zb < -/lmin,glide & Vbx > K,min) 

then glide (a = — 1) by the control law in Sec. IV. A. 3 (48) 
else flap {a = 1) by the control law in Sec. IV. A. 1 and IV. A. 2. 

The switching logic ensures that we have suflicient altitude and forward velocity to glide, but will interrupt 
the constant ascending flapping flight with periods of gliding. 

IV. B. Simulation Results 

Table 1. Key simulation parameters and nominal values 



m=0.3 kg 


h= 


0.0012*eye(3) kgm^ 


R=0.32, c= 0.15 [m] 


k=60 (flapping) or (gliding) 


6^=20 deg 


X= 


10 (flapping) 


A = 30 (ghding) 


Pi=50, p2=30, p3=15 [deg] 


ai = a5=0 


a2 = 


= a6=0 deg 


as = aj = —5 deg 


wing: Cl^Cd'- Eq. (30) 


wing: Cmo = -0.2 


win, 


g: Cma = -0.12 


body: CmO =0.1 


body: Cma = -0.2 



We assume each wing to be a single rigid body. That is, in Eq. (26), we set {xw{r),yw{r), Zw{r))^ = 0. 
Also, we assume that the aerodynamic coefficients do not vary along the chord of the wing blade element. 
This dynamic model with three dimensional wing kinematics is constructed in Matlab/Simulink, allowing 
us to demonstrate how stability can be obtained for flapping flight driven by a CPG network (see Fig. 1). 
A simulation result of alternating flight modes by the control laws in Sec. IV. A is shown in Figs. 7 and 
8. The control parameters and the wing angles driven by the CPG network are shown in Fig. 8. The key 
simulation parameters are listed in Table 1. As shown in Fig. 7, the vehicle begins in a gliding mode with 
cr = — 1, transitions smoothly to ascending flapping flight (cr = 1) at t=7.78 s, and then executes a turn at 
t=ll s while stabilizing both the longitudinal and the lateral modes. The vehicle resumes forward ascending 
flapping flight by driving the bank angle to zero at t=16 s, and the ascending flapping flight switches to 
a glide mode as it reaches some desired maximum altitude t=19.1 s (see Zg in Fig. 7a). Remarkably, the 
flapping flying vehicle makes smooth transitions between various flight modes. This should be a testament 
to the potential of our proposed CPG-based control scheme. 

To turn, at the 11 second mark, we set a desired bank angle 06=40 deg for the symmetry breaking 
control law in Eq. (46), accompanied by setting a constant frequency (a; = 0). What is interesting here 
is that the vehicle settles into a nice banked turn at the same flight path angle without correcting the 
frequency of flapping. This is because the change of ps and pr causes the aerodynamic forces to increase in 
spite of a constant co. The bank angle, rate of turn, and qualitative characteristic (e.g., amplitude of body 
pitch oscillation) of the turn are interestingly linked to the scale and shift of the lead-lag coupling (A32) 
and amplitude, but an exact correlation depends on physical and aerodynamical parameters. This must be 
further investigated and understood in order to implement a better nonlinear control law. The roll angle 
through the turn is about 40 deg and the average global yaw rate is about 45 deg/s with oscillations between 
5 and 75 deg/s. This is feasible in light of Hedrick's experimentation.^-'^ The local angle of attack of each 
wing varies along the wing span (r G [0,i^]) and nicely bounded between ± 50 deg at the steady-state as 
seen in Fig. 9. While we have almost all positive angles of attack in the downstroke, creating high forces, 
the upstroke has a nice balance of positive and negative angles of attack. 

Figure 8b shows the resulting oscillatory behavior of the flapping {(j)w)^ pitch (Ow)^ and lead-lag motion 
(ipw) commanded by the CPG network and highlights the effects of our changing control variables on CPG 
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behavior. From arbitrary initial conditions, the CPG network synchronizes globahy and exponentiahy. 
During ghding flight, the synchronization between the wing joints is less critical, so we set the coupling gain 
A: = 0. Even if the phase difference Aij and the amplitude oscillation pi are time- varying, the synchronization 
errors exponentially tend to zero due to Corollary 2 and Theorem 2. This simulation results supports our 
claim that our CPG-based flapping flight control can stabilize both the longitudinal and lateral dynamics 
by commanding three control parameters: c<;(t), A32 = A76, and ps 7^ pr- 

V. Conclusion 

We investigated the hypothesis that the phase control and synchronization of coupled nonlinear oscillators, 
inspired by central pattern generators (CPGs) found in animal spinal cords, can effectively produce and 
control stable flapping flight patterns and can be used to stabilize the flapping flying vehicle dynamics. An 
engineered CPG network, which ensures the stability and robust adaptation of motion, can signiflcantly 
reduce the complexity associated with engineered flapping flight. Central to the agile flight of natural 
flyers is the ability to execute complex synchronized three-dimensional motions of the wings. In this paper, 
we introduced a mathematical and control-theoretic framework of CPG control theory that enabled such 
synchronized wing maneuvers. In contrast with an independent control of each wing joint by a sinusoidal 
function, the proposed CPG-based method used the phase differences between multiple wing oscillators as 
control mechanisms. We also showed that the central ffight controller, similar to the brain of an animal, can 
stabilize the full 6-DOF vehicle dynamics by commanding a reduced number of control variables such as the 
frequency and phase difference of the oscillators instead of directly controlling multiple joints. Furthermore, 
to our knowledge, this paper presented the first result on alternating fiight modes of gliding and fiapping 
fiight by using the Hopf bifurcation. 

In order to show the effectiveness of the proposed CPG-based flapping flight control, we presented numer- 
ical simulation results by using a realistic vehicle dynamic model with three dimensional wing kinematics. 
This dynamic model includes a tilted stroke plane angle, the lead-lag motion, and the flapping and pitch 
angles of each wing, along with the full states of 6D0F nonlinear vehicle dynamics. We also showed that 
CPG-based flight control can stabilize and alternate two different flight modes of flapping and gliding flight 
by using the synchronized and symmetry-breaking oscillatory motions of two main wings. This result is 
interesting in the sense that the tailless flapping flight dynamics could be effectively controlled without using 
traditional aerodynamic control surfaces. This result agrees with a prior claim of biologists that birds act 
more like tailless aircraft. 

While we understand the challenges associated with lightweight low-power actuators to fully realize the 
potential of three dimensional wing movements, the research described in this paper can further enhance our 
understanding of key mechanisms of biological flyers. Ongoing work includes a study of aeroelasticity with 
various kinds of waveforms and compliant wings. Also, it would be important to identify which wing joint 
variables can be controlled passively to further reduce the complexity. 
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(a) Body velocity and inertial position. 
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(b) Wing joint angles, sync errors, and flapping frequency 
Figure 8. Control inputs and wing joint angles of the two alternating flight modes, flapping and gliding. 
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Figure 9. Local angles of attack of the right wing along the wing span r G [0, R] m with dr = 0.01 m. 
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